Attempted density blowup in a freely cooling dilute granular gas: hydrodynamics 

versus molecular dynamics 



OO 
O 

o 

(N 
C 

Hi 

in 



Andrea Puglisi 1 , Michael Assaf 2 , Itzhak Fouxon 2 , and Baruch Meerson 2 

1 Dipartimento di Fisica, Universita La Sapienza, p.le Aldo Moro 2, 00185 Roma, Italy 

2 Racah Institute of Physics, Hebrew University of Jerusalem, Jerusalem 91904, Israel 

It has been recently shown (Fouxon et al. 2007) that, in the framework of ideal granular hydro- 
dynamics (IGHD), an initially smooth hydrodynamic flow of a granular gas can produce an infinite 
gas density in a finite time. Exact solutions that exhibit this property have been derived. Close to 
the singularity, the granular gas pressure is finite and almost constant. This work reports molecular 
dynamics (MD) simulations of a freely cooling gas of nearly elastically colliding hard disks, aimed 
at identifying the "attempted" density blowup regime. The initial conditions of the simulated flow 
mimic those of one particular solution of the IGHD equations that exhibits the density blowup. We 
measure the hydrodynamic fields in the MD simulations and compare them with predictions from 
the ideal theory. We find a remarkable quantitative agreement between the two over an extended 
time interval, proving the existence of the attempted blowup regime. As the attempted singularity 
is approached, the hydrodynamic fields, as observed in the MD simulations, deviate from the predic- 
tions of the ideal solution. To investigate the mechanism of breakdown of the ideal theory near the 
singularity, we extend the hydrodynamic theory by accounting separately for the gradient-dependent 
transport and for finite density corrections. 
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PACS numbers: 45.70.Qj 



I. INTRODUCTION 



Spontaneous clustering of particles in granular gases 
has attracted much recent interest [l|, H, y, 0, H, [3, 0, H, 
H, E3, EH, E2] • As other pattern-forming instabilities, the 
clustering instability of a freely cooling granular gas has 
served as a sensitive probe of theoretical modeling and, 
first of all, of the Navier-Stokes granular hydrodynam- 
ics (GHD). Although the formal criteria of its validity 
may be quite restrictive, see below, GHD has a great 
power, sometimes far beyond its formal validity limits 
[l3| , in predicting a host of collective phenomena in gran- 
ular flows, such as shocks, vortices and clusters. In the 
recent years, GHD has been applied to a variety of non- 
stationary dilute granular flows [13, El E! EE ES El • 
Non-stationary flows are both appealing and challeng- 
ing for continuum modeling of granular dynamics. As in 
other areas of continuum modeling, this is especially true 
when a non-stationary flow develops a finite-time singu- 
larity E3|. Examples are provided by the finite-time 
blowup of the gas density: at zero gravity [13, EH E3 
(as described by the ideal GHD discussed below, IGHD 
from now on), and at finite gravity [l5[ as described by 
the more complete, non-ideal GHD. Of course, a den- 
sity blowup in a gas with finite-size particles can only be 
an intermediate asymptotics, as the blowup is ultimately 



arrested: either by close-packing effects 11], or by the 
gradient-dependent transport [H. Still, the attempted 
blowup regimes, signaling the development of high den- 
sity regions in the gas, are fascinating and worth a de- 
tailed investigation. One such re gim e has been recently 
addressed by Fouxon et al. [lfA Il7|. They dealt with 
a macroscopically one-dimensional, freely cooling, dilute 
granular flow in the framework of ideal GHD that ne- 
glects the gradient-dependent transport effects: the heat 



diffusion and viscosity. Fouxon et al. derived a class of 
exact solutions to the ideal equations, for which an ini- 
tially smooth flow develops a finite-time density blowup. 
Close to the blowup time r, the maximum gas density 
exhibits a power law behavior ~ (r — i)~ 2 . The velocity 
gradient blows up as ~ — (r — t) -1 , whereas the veloc- 
ity itself remains continuous and forms a cusp, rather 
than a shock discontinuity, at the singularity. The gas 
temperature vanishes at the singularity, but the pressure 
remains finite and almost constant. Extensive numeri- 
cal simulations with the ideal hydrodynamic equations 
showed that the singularity, exhibited by the exact solu- 
tions, is universal, as it develops for quite general initial 
conditions [13, E3| ■ The reason behind this universality 
is in the fact that the sound travel time through the re- 
gion of the developing singularity is much shorter than 
the characteristic inelastic cooling time of the gas in that 
region. As a result, the pressure gradient (almost) van- 
ishes in the singularity region, and the local features of 
this isobaric singularity become essentially independent 
of how the flow was initiated and how it behaves at large 
distances from the singularity. This singularity is of the 
same type as the one that develops, in the framework of 
the IGHD equations, in a general low Mach number flow 
of a freely cooling granular gas [11 . 

Here we perform molecular dynamics (MD) simula- 
tions of a freely cooling gas of nearly elastically colliding 
hard disks, aimed at identifying the "attempted" density 
blowup regime, predicted by the ideal analytical solu- 
tions |17| . We simulate a freely evolving dilute gas of 
nearly elastically colliding hard disks in a narrow channel 
with perfectly elastic side walls. In this geometry both 
the clustering mode in the transverse directions, and the 
shear mode are suppressed (see Refs. @, [3, E3, E3 for 
detailed criteria). As a result, the coarse-grained, or hy- 
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drodynamic flow depends only on the longitudinal coor- 
dinate al ong the channel (and time), as it was assumed 
in Refs. [id, Il7j|. We choose the initial conditions of the 
MD simulations so that the coarse-grained density, ve- 
locity and temperature fields are those producing one of 
the exact blowup solutions of the ideal GHD equations. 
Then we follow the time history of the hydrodynamic 
fields in the MD simulations and compare it with that 
predicted by the ideal exact solution and with numerical 
solutions of non-ideal GHD equations. 

The remainder of this paper is organized as follows. 
In Section II we briefly summarize the ideal GHD anal- 
ysis [ill H3] of the density blowup: we present the ideal 
GHD equations, one of their exact solutions, its main 
features and expected limits of its validity. In Section 
III we describe our MD simulations and compare the hy- 
drodynamic quantities, computed from the simulations, 
with the exact solution of the ideal GHD equations. We 
find that the exact solution is in remarkable quantitative 
agreement with the MD simulations over an extended 
time interval, proving the existence of the attempted den- 
sity blowup regime. As the attempted singularity is ap- 
proached, the hydrodynamic fields, as observed in the 
MD simulations, deviate from the predictions of the ex- 
act solution. To investigate the mechanism of breakdown 
of the ideal solution, we extend the hydrodynamic the- 
ory, in Section IV, in two separate ways. In the first one 
we take into account the gradient-dependent transport: 
the heat diffusion and viscosity, but continue to assume 
that the gas is dilute. In the second one we neglect the 
gradient-dependent transport but take into account, in 
a semi-phenomenological way, finite density corrections. 
Section V summarizes our results and puts them into a 
more general context of hydrodynamic scenarios of clus- 
tering in freely evolving granular gases. 



II. HYDRODYNAMIC THEORY AND DENSITY 
BLOWUP 



Ideal granular hydrodynamics and exact 
solution 



We consider a two-dimensional granular gas of identi- 
cal hard and smooth disks with diameter a and mass set 
to unity, and adopt a simple model where the inelastic 
particle collisions are characterized by a constant coeffi- 
cient of normal restitution r. Throughout this paper we 
will only deal with nearly elastic collisions, 

l-r<l, (1) 

and assume a very small Knudsen number: 



I free/ L <C 1 



(2) 



Here I free is the mean free path of the particles, and L is 
the characteristic length scale of the hydrodynamic fields 
that may depend on time. In addition, we will assume in 



this Section that the local gas density p is much smaller 
than the close-packing density of disks p c — 2/(v / 3er 2 ): 



pa < 1 



(3) 



The strong inequalities © and §5§ need to be verified a 
posteriori, once the hydrodynamic problem in question 
is solved. The strong inequalities {l}-© enable one to 
employ the well-established equations of Navier-Stokes 
granular hydrodynamics (see, e.g. [H, that deal 
with three coarse-grained fields: the mass density p(x, t), 
the mean flow velocity v(x, t) and the granular tempera- 
ture T(x, t). In a sufficiently narrow channel these fields 
depend only on the longitudinal coordinate x, and the 
hydrodynamic equations become 



dp d(pv) 
dt dx 
' dv dv 



= 0, 



P \-di +V d-x 

dT dT _ 
dt dx 

| voVT /dv 
p \dx 



d{ P T) 
dx 



(4) 



-T 



dx 



ApT 



3/2 



«o d 
p dx 



dT 

dx 



(6) 



where A = ^/tt{1 — ?~ 2 )cL vn = (2^/ira) 1 and kq — 
2/(V7T(7), see e.g. Refs. jflUfl]. Equations (J)-© differ 
from the hydrodynamic equations for a gas of elastically 
colliding disks only by the presence in Eq. of the in- 
elastic loss rate term —ApT 3 / 2 , that describes the pro- 
portionality of the energy loss per particle to the num- 
ber of particle collisions per unit time (proportional to 
pT 1 / 2 ) and to the energy loss per collision (proportional 
to T). This additional term, however, brings a whole new 
physics (and mathematics) into the problem. 

Let us rewrite Eqs. ([I])-© in dimensionless variables. 
We will measure the gas density, te mpera ture and the 
velocity in the units of po, To/2 and -y/To/2, respectively, 
where po and To are some characteristic values of the 
initial density and temperature. The time and distance 
will be measured in the units of 



ApoVTo 



and I 



(7) 



respectively. As one can see from Eq. ©, t is the char- 
acteristic cooling time of the gas due to the collisional 
energy loss, while I is the characteristic distance a sound 
wave travels during time r. The numerical factors in 
Eqs. ([7]) are chosen for convenience. We will keep identi- 
cal notation for the rescaled and physical quantities, and 
take care that no confusion arises. Using the rescaled 
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quantities, we rewrite Eqs. 
d(pv) 



dp 
dv 
dT 

i - 



dx 
dv 
dx 
dT _ 
dx 
r 2 d 



0, 



d{ P T) l-r 2 d 



V2p dx 



dx 4y/2 dx 

- T ^-2V2pT 3 / 2 
dx 

VT— 

dx 



T— 

dx 



(l-r 2 )VT 
4^2p 



(8) 
,(9) 



m -(io) 



dv 
dx 



Let us assume that the characteristic magnitudes of the 
rescaled hydrodynamic fields, and of their spatial and 
temporal derivatives, are of order unity (this assumption 
needs to be checked a posteriori). Then we can neglect 
the viscous and thermal conduction terms, as they scale 



as 1 



<C 1, and arrive at the IGHD equations: 



dp 
dt 

dT 
~dt 



d(pv) 



= 0, 



+ v 



dT dv 

= -T 



dv 



dv 



P[ dt 



dx 
2v / 2pT 3/2 



d(pT) 

dx 



,(n) 



(12) 



dx dx 

These equations were investigated in Refs.[l3,0], where 
a family of exact analytic solutions was derived. Here 
we will consider a representative and simple particular 
solution that evolves from the following initial conditions: 



p(x, t = 0) = cosh 



T(x,t= 0) = 2, 



v(x, t = 0) = — 2 arcsin (tanh x) . 



(13) 
(14) 



To remind the reader, we are using rescaled variables 
here. Back in the physical variables, the initial profiles 
are 



Po 



P(x ' 0) = cosh(*/0' 
T(x,0)=T , 

v(x, 0) = — \J 2Tq arcsin 



tanh 



(t) 



(15) 
(16) 
(17) 



where I is defined in Eq. ([7]). That is, at t = the density 
profile has a maximum po and width /, the temperature 
To is uniform, and there is an inflow of the gas towards 
the origin with v(x — > ±oo, t = 0) = ^ir^T a /2. The 
initial scale of variation of the fields, I, is by a factor 
1/(1 — r 2 ) greater than the mean free path of the gas, 
justifying the use of the hydrodynamic description. Fur- 
thermore, both the magnitudes, and the spatial scales of 
the rescaled fields are of order unity which justifies, at 
least for finite times, the use of the ideal equations (jTTJ) 
and ([12")) . Figure [1] shows the initial density and velocity 
fields of the flow in the Eulerian coordinate. 

Now we go over from the Eulerian coordinate x to the 
Lagrangian mass coordinate m = J Q p(x' ,0) dx' , see e.g. 
[211 ] . For the initial density profile (flB"]) . the Eulerian 
coordinate x is related to the Lagrangian coordinate m 
as follows: 



1 / 1 + sin to 
2(m, t = 0) = - In : 

2 VI — sin 777 



(18) 




FIG. 1: The initial values of the hydrodynamic fields in the 
example of attempted density blowup considered in this work. 
Shown are the rescaled initial density and velocity of the gas 
[see Eqs. (|15[) and (|17[) ] versus the rescaled Eulerian coordi- 
nate x/l along the channel. Only the right half of the system 
is shown. The values of po, To and I, used in our molecu- 
lar dynamic simulations, are presented in the beginning of 
subsection III.D. 



Note that the infinite Eulerian interval — oo < x < +oo 
corresponds, in this example, to a finite interval of m: 
— 7r/2 < ?77 < 7r/2, that is to a finite (rescaled) total mass 
of the gas, equal to 7r. In the Lagrangian coordinates, 
the ideal equations ITT|) and (fT2|) arc 



d_ n 
dt Vp 



dv dv 
dm ' dt 



dp 
drn 



|=-2^-2^p3/ V / 2 
dt am 



(19) 
(20) 



where the dilute gas pressure p = pT has been used in- 
stead of the temperature. The exact solution in the La- 
grangian coordinates is the following [l6|, [ItJ : 



COS 777 



p(m,t) = 

(1 — tcosm) z 

?j(777, t) = — 2777 + 2t Sin 777 



p{m, t) — 2 cos 777, (21) 
(22) 



As x(m,t) can be calculated explicitly (l7| . 



x(m, t) 



In 



1 + sin 777 
1 — sin to 



-2tm + t 2 sin m, (23) 



equations (f2"Tj) - (r2"3"l) describe the time history of the hy- 
drodynamic fields in the x-coordinate in a closed para- 
metric form. Let us consider some important features of 
this simple exact solution. 



B. Density blowup and its properties 

A momentary look at the first of Eqs. ([2T]) reveals that 
the density at the origin blows up in a finite time. Back 
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in the dimensional variables one obtains 

Po 



p(x = 0,t) 



{l-t/rf 



(24) 



The gas temperature becomes zero at the singularity, the 
velocity gradient dv(x,t)/dx blows up as ~ — (r — t) -1 , 
whereas the velocity itself remains continuous and forms 
a cusp. This behavior at singularity is quite different 
from that of the free flow singularity (that one observes 
for the same initial velocity profile but a zero gas pres- 
sure). In particular, for the free flow singularity the den- 
sity blows up as (1 — t/r)~ l , while the velocity develops 
a shock discontinuity [22|]. See Refs. [la, GjJ for more 
differences between the two types of the singularities. 

It was found in Refs. [HI, [TtJ that the finite time 
blowup of the gas density and of the velocity gradient is 
not a consequence of specially chosen initial conditions. 
Rather, it appears, in the framework of the ideal equa- 
tions (fTI) and (fT2"|) [or, equivalently, of Eqs. ([il?]) and 
(|2U)) ]. for quite general initial conditions. A robust lo- 
cal feature of this singularity is a finite, non-zero value 
of the gas pressure at the time of the density blowup. 
The singularity is universal because the sound travel time 
through it is much shorter than the characteristic inelas- 
tic cooling time of the gas. As a result, the pressure 
gradient (almost) vanishes in the singularity region, and 
the local features of this isobaric singularity become in- 
dependent of the flow details at large distances. 

As an infinite density cannot be reached in a gas with 
finite-size particles, it is clear that, sufficiently close to 
the attempted singularity, some of the assumptions made 
on the way to the ideal equations (jUJ)- (H3) break down. 
However, independently of the precise way of regulariz- 
ing the infinite density, an attempted singularity implies 
the formation of a region with a very high density. Fur- 
thermore, the ideal solution should accurately describe 
the evolution of the system for times not too close to the 
attempted density blowup time r. Before we look into 
where the ideal theory breaks down, let us consider some 
global characteristics of the exact solution. For these the 
singularity time t = t turns out not to be special. First, 
we note that the solution describes a gas with a (con- 
stant) finite number N of particles, given by 



N = L, 



p(x)dx — irpoLyl — 



2j2nL 



y 



(1 



(25) 



where L y is the channel width. Note that N is indepen- 
dent of po- for larger po the initial density profile has 
a higher peak but a smaller width, so that iV remains 
constant. Another global characteristics of the solution 
is the total energy of the gas 



E(t) = Ly 



pv 



pT I dx . 



For the thermal part of the energy we find, after a simple 



algebra, 



Etnit) 



pTdx 



2 I pTdx 
Jo 



tt/2 



= 2p T Q l / [1 - (t/r) cos m} 2 dm 
Jo 

2" 

I I \ :-: I I \ 

7T - 4 



PqTqI 



7T / t 

2 [t 



(26) 



Similarly, for the kinetic energy of the macroscopic mo- 
tion we find 



Ekin (t) 



= PqTqI 



12 It/ 2 It 



(27) 



Summing the two and using Eq. (f2"5"l) , we obtain (Tt| 



E(t) = NT Q 



7T 4 8 ft 
12 ~ n It 



(28) 



As time increases, E(t) is monotone decreasing for t < 
t. We also observe that t = r is a regular point of 
E(t), where nothing dramatic happens. Note that the 
parabolic-in-time law of the energy decay is quite differ- 
ent from Haff's law E(t) = E(0) /(l + 2t/r) 2 obtained for 
freely cooling homogeneous granular gas with density po 



C. Breakdown of ideal theory 

For a one-dimensional flow, the applicability of the so- 
lution (T2"Tj) - (f2"5)) is determined by the applicability of the 
ideal equations (|11[) - ([T2"|) that it solves exactly. Analysis 
in Ref. [17J shows that, sufficiently close to the attempted 
singularity, the ideal equations become invalid. This hap- 
pens because of one of two reasons (or both): either the 
gas becomes dense, so that criterion ([3]) breaks down, or 
the heat diffusion becomes important, invalidating the 
ideal equations (flT|) - p^)) . The time % r at which the ideal 
theory breaks down can be estimated as follows 17] : 



1 — — max(V Poo- 2 , 1 — r). 



(29) 



Therefore, the "bottleneck" for the validity of the equa- 
tions is set by the initial conditions: if the maximum 
is determined by the first (correspondingly, the second) 
term in the right hand side of Eq. (|29|) . the ideal equa- 
tions become invalid because of the finite gas density 
(correspondingly, the finite heat diffusion). As each of 
the two terms is very small by assumption, the solution 
is expected to break down only close to the attempted 
singularity (24|. 

The one-dimensional solution may also become invalid 
because of instability with respect to small initial pertur- 
bations that are inevitably present in MD simulations. 
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Numerical solutions of the hydrodynamic equations, re- 
ported in Refs. 0, 0], strongly suggest that the ideal 
solution is stable with respect to small longitudinal per- 
turbations. This does not exclude possible instability 
with respect to small transverse perturbations. The only 
available analytic result here is the one obtained from 
the stability condition for a homogeneous cooling state, 
see Refs. @, H, [1(1 HU- That stability criterion comes 
from a competition between the (destabilizing) inelastic 
cooling and the (stabilizing) heat diffusion and viscosity 
in the transverse direction. The stability criterion de- 
mands that L y be less than a threshold value depending 
on 1 — r, po and a. The stability problem for the strongly 
inhomogeneous and time-dependent exact solution is ob- 
viously more complicated, and its complete analytic so- 
lution does not seem feasible. It is therefore important 
that our MD simulations, presented in the next Section, 
strongly suggest that, for sufficiently narrow channels, 
no instability in the transverse direction occurs for the 
time-dependent flow we are working with. 

Let us summarize the main theoretical predictions. For 
the initial conditions (fl~3| and (fl"4|) . a nonlinear time- 
dependent flow sets in, described by the ideal exact solu- 
tion: Eqs. (f2"T ]) -([2"3"] ) . This flow "attempts" to develop a 
density blowup. However, close to the attempted singu- 
larity one (or both) of the two factors: the finite density 
and the heat diffusion, invalidates the solution. The rel- 
ative importance of the two factors is determined, via 
Eq. ([29]) , by the initial conditions. 



III. MOLECULAR DYNAMICS SIMULATIONS 
A. General 

To test the theoretical predictions, we performed MD 
simulations of a free cooling granular gas in a narrow two- 
dimensional channel. The initial conditions correspond 
to hydrodynamic profiles (fT3|) and (fT4|) and satisfy the 
strong inequalities (H}-©. According to the theory, they 
are expected to generate the nonlinear time-dependent 
flow described by Eqs. ([2l| -(|23 ]) . Our MD simulation 
calculates the evolution of a gas of N' identical inelastic 
hard disks of unit mass, with diameter a, in a channel of 
dimensions L' x = L x /2 and L y . As the expected hydro- 
dynamic flow is symmetric with respect to x = 0, only 
one half of the system, x E [Q,L x /2], is simulated, so 
N' = N/2. Each wall of (the one half of) the channel 
is solid and reflects elastically the disks colliding with it. 
The particles move freely until a collision ("event") oc- 
curs when two disks i and j find themselves at a distance 
equal to a. The collision is resolved instantaneously, leav- 
ing the positions of the particles unaltered and updat- 
ing their velocities from (vj,Vj), before the collision, to 
(v-,Vj), after the collision. The update rule conserves 
the total momentum and reduces the total kinetic energy, 



with a constant coefficient of normal restitution r € [0, 1]: 

v' t =v 4 -i±I(g.<x)<x, (30) 
v' j = v j + ±±l(j B >&)&, (31) 

where g = v< — Vj and <x is the unit vector joining the 
centers of the two disks. The hard-core interactions make 
possible the following optimization of the algorithm. It 
is sufficient to calculate the first collision times of all 
particles and then select the absolute first one. The 
system is freely evolved up to that time, then the col- 
lision is resolved, and a new list of collision times is 
computed. With standard optimization techniques of the 
search procedure it is possible to achieve fast computa- 
tion times [25]. Nevertheless, the time performance is 
proportional to the number of collisions occurred, so the 
ratio between the physical time and the cpu-time goes 
down when dense clusters emerge in the system. 



B. Initial conditions 

The initial position and velocity of each of the N' disks 
are chosen randomly with probability distributions corre- 
sponding to the desired initial hydrodynamic fields. This 
was implemented with the following procedure. For each 
disk i 

1. the longitudinal position Xi is chosen with probabil- 
ity proportional to p(x, 0) from Eq. ([T5]) for x > 0, 
using an acceptance/rejection method: 

(a) a random ^-position is generated with uni- 
form probability on the interval [a/2,L' x — 

(b) a random number z with uniform probability 
on [0, m.ax{p(x, 0)}] is compared with p(xi, 0), 

(c) if z < p(xi,0) the position is accepted, other- 
wise the procedure is repeated from Hal 

2. then the vertical position yi is chosen with uniform 
probability on the interval [o~/2, L y — cr/2], 

3. a non-overlap check is performed: the distance be- 
tween the disk center (xj, yi) and all the previously 
placed disk centers must be greater then a: if the 
condition is not satisfied, the procedure is repeated 
from Hal 

4. the velocity components vf and vf are chosen from 
a Gaussian distribution with zero mean and vari- 
ance equal to To; then the longitudinal component 
is shifted by an amount v(x, 0) from Eq. (|17p with 
x > 0. 



G 



C. Lagrangian coordinate and hydrodynamic fields 

We verified that, for our choice of the channel dimen- 
sions, the gas remained homogeneous in the y-direction. 
By virtue of this observation, it was sufficient to deal with 
one-dimensional hydrodynamic profiles depending on x. 
For a direct comparison with the analytical solution of 
the IGHD, the hydrodynamic profiles were obtained us- 
ing a uniform binning in the Lagrangian mass coordinate. 
Using the same notation as in Section I, we define the La- 
grangian mass interval for the simulated flow as [0, 7r/2], 
where 7r/2 corresponds to (one half of) the total gas mass, 
N 1 = N/2. Let nbi n be the number of bins chosen to sam- 
ple the hydrodynamic profiles and N b in — N' /nun be the 
average number of particles per bin. At a given time t 
all particles are ordered so that Xi < £j+i, i € [0, N' — 1]. 
Then each bin j G [l,7i{,j n ] has its leftmost border at 
x [{j-i)N bin } an d its rightmost border at £[jAr bi „-i] ■ These 
bins are non-uniform in the x coordinate, but are uni- 
form in the Lagrangian mass coordinate as each contains 
the same mass N b i n . The position of the j-th bin is 



m ; 



*ti ~ V2) 



All particles belonging to the j-th bin contribute to the 
value of the hydrodynamic fields: 



p(mj,t) 
v(mj,t) 

%,<)= Ei£j 



LyLj 

Nbin 



2N br , 



(32) 
(33) 

v 2 (m j: t),(34) 



where we have used the shorthand notation i € j to 
denote particles in the j-th bin, and Lj to denote the 
length of the j-th bin. The pressure field p(m,-,i) = 
p(rrij,t)T(mj,t) is obtained straightforwardly. 

The hydrodynamic fields, computed for individual re- 
alizations, exhibit a significant noise. To get rid of the 
noise, all hydrodynamic profiles presented in the next 
Section were obtained after an averaging over fOO MD 
simulations with different initial conditions, correspond- 
ing to the same initial hydrodynamic fields and obtained 
with the procedure described in subsection B. 



D. MD simulations versus ideal solution 

The following parameters were chosen for the simula- 
tions: a = 1, po = 1CT 4 , T = 1, N' = 5 x 10 4 , and 
L y = 125. For convenience, the coefficient of normal 
restitution r was chosen so that 1 — r 2 = \/tt/2 x 10~ 2 , 
i.e. r — 0.99371367 .... This choice of parameters sets 
I ~ 1.27324 x 10 6 , L' x = 10/ and r ~ 1.8006 x 10 6 . The 
evolution of the density field, as obtained in the simula- 
tions and as predicted by the ideal theory, is displayed 
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FIG. 2: (Color online) Evolution of the density p(m,t) in 
the Lagrangian frame at initial times. The circles: the re- 
sults of MD simulations. The solid line: the prediction of 
the analytic solution, first of Eqs. (|21[l . back in the physi- 
cal variables. The dashed line: a numerical solution of the 
non-ideal hydrodynamic equations (f8)l- (|10[) that account for 
the gradient-dependent transport. The dashed line is indis- 
tinguishable from the solid line. 



in Fig. [5] for times up to t — 0.8t, and in Fig. [3] for 

later times, up to time t — 1.055r. The figures show 
that the ideal solution is in remarkable agreement with 
the MD simulations up to times t ~ 0.9r. At later times, 
when the density peak exceeds ~ 10 -2 , the ideal solution 
starts to deviate from the MD simulation in the neigh- 
borhood of m = 0. The actual density peak continues 
to grow, but slower than predicted by the ideal solution. 
At time t = r, when the ideal solution predicts the den- 
sity blowup in x = 0, the actual density p(0,t) ~ 0.2. 
The close-packing density p c = 2/(\/3cr 2 ) is reached at 
t ~ 1.05r, see the last frame of Figure[31 Sufficiently far 
from m = 0, the ideal solution remains very accurate. 
(We checked that this statement remains true even be- 
yond the attempted singularity time: until the end of the 
MD simulations.) 

The gas velocity profiles, shown in Figs. [5] and O in 
a linear and logarithmic scale, respectively, are very ac- 
curately predicted by the ideal solution, Eq. (|2"T]) . until 
late times. Surprisingly, the excellent agreement remains 
even at times greater than 0.9t, when the density peak 
already significantly deviates from the theoretical one. 
To be able to see the small deviations from the theory, 
we had to use, in Fig. O a logarithmic scale. 

Similarly, an inspection of the pressure profiles, see 
Fig. [6l shows an excellent agreement with the prediction 
of the ideal theory, p(m, t) = PqTq cos(m), see the second 
of Eqs. (|2"2"|) . Discrepancies of about 2% appear only at 
late times, when the density is already about one half of 
the theoretically predicted value. At times close to r, the 
pressure field, as found in the MD simulations, develops 
a dip close to x = m = 0, reaching a value about 15% 
lower than the theoretically expected value p(0) = 10~ 4 . 
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FIG. 3: (Color online) Evolution of the density p(m, i) in the Lagrangian frame at late times. The m-axis is zoomed in, 
in order to focus on the high density region. The circles: the results of MD simulations. The solid line: the prediction of 
the first of Eqs. (|21|l . back in the physical variables. The dashed line: a numerical solution of the non-ideal hydrodynamic 
equations (|8])- (|10[) that account for the gradient-dependent transport. The dash-dotted line marks the close packing density 
p c = 2/(V3a 2 ). 
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FIG. 4: (Color online) Evolution of the longitudinal velocity 
v(m,t) in the Lagrangian frame. The circles: the results of 
MD simulations. The solid line: the prediction of Eq. (122)1 . 
back in the physical variables. The dashed line: a numer- 
ical solution of the non-ideal hydrodynamic equations ([8}- 
(|10|) that account for the gradient-dependent transport. The 
dashed line is indistinguishable from the solid line. 



FIG. 5: (Color online) Evolution of the (minus) longitudinal 
velocity v(m, t) in the Lagrangian frame, in a semi-logarithmic 
scale. The circles: the results of MD simulations. The solid 
line: the prediction of Eq. (I22|l . back in the physical vari- 
ables. The dashed line: a numerical solution of the non- 
ideal hydrodynamic equations ([8)l- (|10[l that account for the 
gradient-dependent transport. 



A direct characterization of the attempted gas density 
blowup is provided by the time history of the density at 
x — 0. The ideal solution predicts, see Eq. (|24j) . that 

i_ n*r=*. ( 35) 

This prediction is extremely well supported by the MD 



simulations in Fig. [71 until t ~ 0.9r. The subsequent 
deviation from the theory appears as saturation of the 
quantity 1 — y/ po/ p(0,t) at the value 1 — y/po/pc cor- 
responding to the close packing density p c . The same 
Fig. [7] also depicts a different quantity, 1 — \/T(0, t). In 
view of the theoretical expectation p(0, t) = pqTq = p$, 
this quantity is also expected to grow as t/r, and Fig. [7] 
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FIG. 6: (Color online) Evolution of the gas pressure, obtained 
using the ideal equation of state p(m,t) — p(m,t)T(m,t), in 
the Lagrangian frame. The circles: the results of MD simula- 
tions. The solid line: the prediction of the second of Eqs. (121[) . 
back in the physical variables. The dashed line: a numeri- 
cal solution of the non-ideal hydrodynamic equations (l8|)- (110[) 
that account for the gradient-dependent transport. 




FIG. 7: (Color online) The time history of the density maxi- 
mum p(0, i) and the temperature minimum T(0, i). The sym- 
bols show the results of MD simulations: black circles depict 
the quantity 1 — *J po/p(0, t), green diamonds depict the quan- 
tity 1 - y/T(0,t). The solid line is the prediction of Eq. (gS) : 
an immediate consequence of Eq. (|24jl . The dashed line is a 
numerical solution of the non-ideal hydrodynamic equations 
((H))- (|lUp that account for the gradient-dependent transport. 
The inset zooms in at later times. 



indeed shows this growth until t ~ 0.9r. The quantity 
1 — -\/T(0, t) also saturates at later times, but at a value 
slightly different from 1 — \J po/p c , see the inset of Fig. [7] 
This is consistent with the pressure deviation from po at 
very late times. 

We also present, in Fig. [51 the time history of the total 
energy per particle, 



N' 



N' 



N' ^ 2 



(36) 




FIG. 8: (Color online) Decay of the total energy per particle 
E(t)/N' as measured in the MD simulations (the circles) and 
predicted by the exact solution (the solid line). The time 
derivative (N')~ 1 dE/dt, displayed in the inset (the notation 
for the circles and line is the same), shows that the energy 
decay remains smooth at all simulated times. 



as found in the MD simulations, and compare it with the 



theoretical prediction, Eq. (|28j) . Here the agreement is 
very good at all times, with a 3% error at late times. 
The (numerical) time derivative of E(t), depicted in the 
inset of Fig. [SI remains smooth also close to the singular- 
ity time t. Actually, this is not surprising, as the main 
contribution to the thermal energy of the gas is made by 
the peripheral gas (in the Lagrangian frame), which is 
hotter and more dilute than the gas in the region close 
to ttl = 0. Furthermore, the main contribution to the ki- 
netic energy of macroscopic motion is again made by the 
peripheral gas (in the Lagrangian frame) which moves 
faster than the gas in the region close to m = 0. The 
peripheral gas continues to follow the ideal theory at all 
simulation times, and this explains the remarkable suc- 
cess of the ideal solution in predicting the total energy 
history. 

To conclude this Section, the MD simulations clearly 
show, over an extended period of time, the existence of 
the "attempted" density blowup regime. The ideal gran- 
ular hydrodynamics (IGHD) predicts very accurately the 
hydrodynamic profiles, observed in the MD simulations, 
up to times close to the attempted singularity. The den- 
sity field, as measured in the MD simulations, starts to 
deviate from the ideal theory at time t ~ 0.9r. Some- 
what surprisingly, the rest of the hydrodynamic fields 
continue to show good agreement with the theory un- 
til even closer to the attempted singularity time r. In 
the following Section we will see that the agreement with 
theory at later times improves significantly when the non- 
ideal hydrodynamic equations (|5j)- (fTT))) . that account for 
the gradient-dependent transport, are employed. 



IV. NON-IDEAL HYDRODYNAMICS 

To investigate the mechanism of breakdown of the ideal 
theory, we extended the hydrodynamic theory in two 
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become significant, and the NIGHD profiles approximate 
the MD simulation results much better than the ideal the- 
ory. As the maximum density continues to increase (and 
finally approaches the close packing density p c ), the di- 
lute NIGHD description ultimately breaks down. In Fig. 
[3] it occurs at about t/r ~ 0.98. 

In the second type of hydrodynamic computations we 
discarded the viscous and heat diffusion terms but took 
into account (moderate) finite density effects. This was 
done by adopting, in Eqs. (fl"§|) and (|2"0)) . instead of the 
ideal equation of state and ideal energy loss rate, the 
Carnahan-Starling equation of state [27] and a modifi- 
cation of the energy loss rate, derived by Jenkins and 
Richman [28| in the spirit of Enskog theory: 




FIG. 9: (Color online) The density profiles for t/r = 0.944 
(a) and t/r = 0.966 (b) as obtained from the MD simulations 
(the circles), the analytic solution (the dashed-dotted line), 
the numerical solution of the dilute NIGHD equations (the 
dashed line), and the numerical solution of finite-density hy- 
drodynamic equations with the gradient-dependent transport 
terms neglected (the solid line). 



separate ways. In the first one we took into account 
the gradient-dependent transport: the heat diffusion and 
viscosity, but continued to assume the the gas is dilute. 
In the second one we neglected the gradient-dependent 
transport but took into account finite density corrections. 
The hydrodynamic equations were solved numerically in 
Lagrangian coordinates, using an accurate variable mesh 
and variable time step solver [26| . 

First, we solved (the Lagrangian form of) Eqs. (|8|)- (fl0|) 
of non-ideal granular hydrodynamics (NIGHD). These 
equations account for the viscous and heat diffusion 
terms but still assume a dilute gas. 

The numerically obtained NIGHD profiles are pre- 
sented, together with the MD simulations and the ideal 
analytical solution, in Figs. [2][7j As expected, the 
NIGHD profiles coincide with the MD simulations and 
with the ideal theory for early and intermediate times. 
At late times, the gradient-dependent transport terms 



P 
A 



PT 



(37) 



where 



g(p) 



i - 



7irp 
32V3 



1 



7rp 
2V3 



is the equilibrium pair correlation function of hard disks 
at contact. In the dilute limit p — > one obtains g = 1 
and recovers the ideal equation of state p = pT . 

In Fig. [5] we compare four different results for the gas 
density: the MD simulations, the ideal analytical solu- 
tion, the numerical solution of the first type (the NIGHD- 
equations), and the numerical solution of the second type. 
It is clearly seen that, for the choice of parameters used in 
our MD simulations, the numerical solution of the second 
type is not as successful as that of the first type. This 
could be expected, as for the times when the maximum 
density is still much smaller than the close packing den- 
sity, the finite-density corrections (which are of the order 
of pj p c ) are still negligible. In contrast, the numerical 
results from the NIGHD equations agree well with the 
MD simulations and show that, as the attempted density 
blowup is approached, the heat conduction and viscosity 
effects can become important when the gas density is still 
small. 



V. SUMMARY AND DISCUSSION 

Our MD simulations proved the existence of an at- 
tempted density blowup regime as described by an exact 
solution of the ideal granular hydrodynamic equations. 
We found the ideal solution to be in remarkable quan- 
titative agreement with the MD simulations over an ex- 
tended time interval, but not too close to the attempted 
singularity. As the attempted singularity is approached, 
the exact solution breaks down. A more complete hy- 
drodynamic theory, that accounts for the heat diffusion 
and viscosity, but still assumes a dilute gas, continues 
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to agree with the MD simulations until the gas density 
becomes a fraction of the close packing density of disks. 

Let us put the results of this work into a more general 
context of clustering instability of a freely cooling granu- 
lar flow. As we have already noted, the local properties of 
the density blow-up, exhibited by the exact solutions of 
the IGHD equations, are the same as the local properties 
of the density blow-up exhibited by a low Mach number 
flow of a freely cooling granular gas with the heat dif- 
fusion neglected [Hj]. A low Mach number flow emerges 
when the pressure balance sets in on a shorter time scale 
than the temperature balance. In this case any local in- 
elastic cooling causes a (low Mach number) gas inflow 
into the colder region so as to increase the local gas den- 
sity there and keep the pressure gradient (almost) zero. 
The resulting density instability develops on the back- 
ground of an (almost) homogeneous gas pressure. This 
is consistent with the MD simulations presented here, see 
Fig. [6l where the pressure in the vicinity of the density 
maximum hardly changes up to times very close to the at- 
tempted singularity time. For brevity we will call the low 
Mach number flow instability scenario Scenario 1. Sce- 
nario 1 first appeared in astrophysics and plasma physics 
in the context of condensation instabilities in gases and 
plasmas that cool by their own radiation [29(. 

As many as four additional hydrodynamic scenarios 
of clustering in a freely cooling granular gas have been 
discussed in the literature. We start with the pressure 
instability scenario, or Scenario 2. It was discussed, in 
the context of the granular clustering, by Goldhirsch and 
Zanetti Q , although it was also known earlier to the as- 
trophysics and plasma physics communities, see Ref. [29| 
for a review. Scenario 2 is usually presented in the fol- 
lowing way. Let us consider a small local increase in the 
gas density. This increase causes an increase in the col- 
lisional energy loss. As a result, the gas pressure falls 
down, a gas inflow develops, causing a further density 
increase, and the process continues. Importantly, Sce- 
nario 2 assumes that the inelastic cooling time is much 
shorter than the sound travel time. In other words, it is 
the local temperature balance that sets in rapidly here, 
and the resulting pressure gradient drives the flow on a 
relatively slow time scale. 

As of present, there has been no detailed nonlinear 
analysis behind Scenario 2. The (well established) linear 
stability theory of the homogeneous cooling state [H, Q 
indicates that Scenarios 1 and 2 operate in two oppo- 
site limits: for sufficiently short and long perturbation 
wavelengths, respectively. The physics behind this is the 
following. The characteristic cooling time due to the in- 
elastic collisions is independent of the length scale of the 
initial perturbation, whereas the sound travel time scale 
is proportional to it. As a result, when all other parame- 
ters are fixed, Scenario 1 corresponds to an intermediate 
wavelength limit of the clustering instability, while Sce- 
nario 2 corresponds to the long wavelength limit. (In the 
short wavelength limit the homogeneous cooling state of 
the gas is stable, as the clustering instability is suppressed 



by the heat diffusion 

Now let us consider Scenario 3 that also assumes a 
long- wavelength limit. As the gas temperature falls down 
rapidly because of the inelastic cooling, the flow is de- 
scribable by the zero pressure (or flow by inertia) ap- 
proximation [l(| ■ Were it not regularized by close pack- 
ing effects, such a flow would develop a finite-time density 
blow up ( of a different type than the low Mach number 
flow) [T(J, [22J ■ If the compressional heating interferes ear- 
lier than the close packing effects, the pressure becomes 
relevant again, and Scenario 3 gives way to the Scenario 
1 [l6|. In the opposite case the late time dynamics of 
the system is describable by the Burgers equation [ill ]. 
Which of the two regimes is realized in a particular set- 
ting depends on the initial conditions. 

Scenarios 1-3 do not invoke the shear mode instability, 
and so they can operate both in one-dimensional, and 
multi-dimensional settings. On the contrary, Scenarios 4 
and 5 do invoke the shear mode, and so they are intrin- 
sically multi-dimensional (and intrinsically non-linear). 
Scenario 4 exploits the fact that the unstable shear mode 
may contribute, via a non-linear coupling, to the growth 
of the clustering mode. Obviously, the nonlinear coupling 
is the only hydrodynamic mechanism of driving the clus- 
tering mode if the system size is larger than the critical 
size for the shear mode instability, but smaller than the 
critical size for the clustering mode instability. Further- 
more, numerical analysis, performed in Ref. [6], indi- 
cated that the nonlinear coupling plays a dominant role 
in the initial density growth also in the case when the 
system size is comparable to the critical system sizes for 
the clustering and shear instabilities. (More precisely, 
the wave number of the monochromatic test perturba- 
tion of the transverse velocity in Ref. 6] was within 
the instability regions of both the shear, and the clus- 
tering modes. However, twice the wave number already 
came out of the instability region.) What happens in 
much larger systems is presently under investigation. It 
turns out that well above the clustering mode instabil- 
ity threshold the nonlinear coupling with the shear mode 
does not dominate the density growth (though it does 
make the theory more cumbersome). In the channel ge- 
ometry, that we adopted in this paper and in the previous 
works [l(| [0, d; [I3i [HI j the shear mode is suppressed, 
and the clustering instability develops in its pure and 
simplest form. 

Now we proceed to Scenario 5 [2| that exploits the 
fact that the unstable shear mode heats the gas in some 
regions. Scenario 5 assumes that this heating can be 
balanced by the inelastic energy loss, rendering (quite 
a complicated) steady state. It is furthermore assumed 
that this steady state is unstable with respect to small 
perturbations, and it is this instability that causes the 
granular clustering. We are unaware of a quantitative 
theory that would support Scenario 5, or of any quantita- 
tive test of Scenario 5 in MD simulations or in numerical 
solutions of hydrodynamic equations. 

To complete the comparison of the five hydrodynamic 
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scenarios of clustering we note that the only scenarios 
that have addressed, up to date, a strongly nonlinear 
stage of the clustering process quantitatively is the low 
Mach number flow instability scenario (Scenario 1) (l8j 
and the zero pressure scenario (Scenario 3) [l(J El- ^ 
is the consideration of a strongly nonlinear stage that 
enables one to identify attempted finite-time density 
blowups: prototypes of the dense granular clusters. 

Which results of this work will withstand a general- 
ization to more realistic granular flow conditions: for ex- 
ample, rotational degrees of freedom and tangential in- 
elasticity of collisions? Including the rotational degrees 
of freedom and tangential inelasticity of collisions into 
a hydrodynamic description is possible under some lim- 
itations [HI, H3|- Solving the corresponding nonlinear 
hydrodynamic equations analytically will of course be 
beyond our reach. It is likely that, when the gradient- 
dependent transport is negligible, these nonlinear equa- 
tions will again exhibit a finite-time density blowup. In- 
deed, the development of closely packed granular clus- 
ters in a granular flow is a robust phenomenon. There- 
fore, it is natural to conjecture, based on results of this 
work, that more realistic granular clusters (those emerg- 



ing when the rotational degrees of freedom are taken into 
account) will still be describable as regularized attempted 
density blowups. 

In summary, the results of the present work gives sup- 
port to the notion of a granular cluster as a regularized 
density blowup of ideal granular hydrodynamic equa- 
tions, put forward in Refs. 0, EE E EE EE Cl- 
in more general terms, they present additional evidence 
that granular hydrodynamics is a powerful and accurate 
quantitative theory of granular flows, especially once it 
is employed within its limits of applicability but, luckily, 
sometimes even beyond them. 
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